Edge states in graphene quantum dots: Fractional quantum Hall effect analogies and 

differences at zero magnetic field 
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We investigate the way that the degenerate manifold of midgap edge states in quasicircular 
graphene quantum dots with zig-zag boundaries supports, under free-magnetic-field conditions, 
strongly correlated many-body behavior analogous to the fractional quantum Hall effect (FQHE), 
familiar from the case of semiconductor heterostructures in high magnetic fields. Systematic exact- 
diagonalization (EXD) numerical studies are presented for the first time for 5 < N < 8 fully 
spin-polarized electrons and for total angular momenta in the range of N(N — l)/2 < L < 150. We 
present a derivation of a rotating-electron-molecule (REM) type wave function based on the method- 
ology introduced earlier [C. Yannouleas and U. Landman, Phys. Rev. B 66, 115315 (2002)] in the 
context of the FQHE in two-dimensional semiconductor quantum dots. Lhe EXD wave functions 
are compared with FQHE trial functions of the Laughlin and the derived REM types. It is found 
that a variational extension of the REM offers a better description for all fractional fillings compared 
with that of the Laughlin functions (including total energies and overlaps), a fact that reflects the 
strong azimuthal localization of the edge electrons. In contrast with the multiring arrangements 
of electrons in circular semiconductor quantum dots, the graphene REMs exhibit in all instances 
a single (0, N) polygonal-ring molecular (crystalline) structure, with all the electrons localized on 
the edge. Disruptions in the zig-zag boundary condition along the circular edge act effectively as 
impurities that pin the electron molecule, yielding single-particle densities with broken rotational 
symmetry that portray directly the azimuthal localization of the edge electrons. 

PACS numbers: 73.21. La, 73.43.Cd 



I. INTRODUCTION 



Since the discovery 1 of the fractional quantum Hall 
effect (FQHE) in two-dimensional (2D) semiconductor 
heterostructures in the presence of a high perpendicular 
magnetic field (B), phenomena associated with strongly 
correlated electrons in the lowest Landau level (LLL) 
have attracted significant and continuous attention. 2-13 
Early on, it was realized that the essential many-body 
physics in the LLL could be most effectively grasped 
through the use of trial wave functions, with celebrated 
examples being the Jastrow-type Laughlin 2 (JL) and 
composite fermion 5 (CF) trial functions associated with 
the formation of a special class of quantum- liquid states. 6 
Later interest in finite 2D electronic systems, like semi- 
conductor quantum dots (QDs) under high B, led to the 
consideration of a different class of analytic trial functions 
known as rotating electron (or Wigner) molecules 8-11 
(REMs or RWMs). An advantage of the REMs is that, 
while they exhibit good total angular momenta, they di- 
rectly incorporate the molecular (crystalline) configura- 
tions that dominate the anisotropic pair correlation func- 
tions revealed through numerical exact-diagonalization 
studies for a finite number of electrons under high B in a 
disk geometry. The initial derivation 8 of the REM trial 
functions generated a flare of theoretical activity around 
the question which class of trial functions (or combina- 
tion of them) is most appropriate for describing the corre- 
lated many-body physics in the LLL of a small number of 
electrons N. 7 > 9 ~ 13 Furthermore, experimental advances in 
the field of ultracold trapped neutral atoms have been fol- 



lowed by considerable theoretical activity regarding the 
nature of correlated states in the LLL that are formed 
during the rapid rotation of the trap; see, e.g., Refs. 14- 
19. 

Recent progress in the fabrication of new materials, 
and in particular in isolating and handling of a single 
graphene sheet, 20-23 offers most promising materials for 
future, post-silicone, miniaturized electronics 24,25 (some- 
time called nanoelectronics) . This expectation is based 
on the two-dimensional character of graphene, where the 
electrons are essentially confined in the spatial direction 
normal to the graphene plane. Fabrication of nanoscale 
device elements for use in electronics, spintronics and 
information processing, such as single-electron transis- 
tors, quantum point contacts, and quantum dots, would 
require additional confinement in the other two spatial 
dimensions. However, to achieve the requested addi- 
tional confinement, techniques (based on electrostatic 
gating) developed for the creation of QDs in semicon- 
ductors (such as GaAs) cannot be used because of the 
unique electronic structure of graphene. The difficulty 
originates from the relativistic, Dirac-like, nature of the 
low-energy quasiparticles in graphene. In particular, the 
gapless nature of the electrons in graphene 26 allows them 
to penetrate unimpeded through a high and wide poten- 
tial barrier. 27 This phenomenon, which is known as the 
Klein paradox, 28 ' 29 is in fact not a paradox but a conse- 
quence of the relativistic character of the electrons, with 
a sufficiently high potential being repulsive for electrons 
but attractive for positrons, thus resulting in positron 
states inside the barrier which can be matched to the 
electronic continuum states outside, consequently result- 
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ing in perfect transmission through the barrier; the un- 
derlying property of the Dirac equation is known as the 
charge-conjugation symmetry. 



In light of the above, one wishes to explore alterna- 
tive, non-electrostatic methods, for fabrication of lower 
dimensionality graphcnc nanostructures. One route for 
achieving the desired added planar confinement consists 
of etching, or cutting graphene into the desired geometry 
(e.g., ribbons, 30-32 circular disks, or other shapes 33-35 ). 
It is expected that further progress in fabrication, char- 
acterization and understanding of the properties of such 
graphene nanostructures (in particular zero-dimensional 
QDs) would lead to their use for the study of interesting 
many-body phenomena, as well as their employment as 
components in miniaturized electronic devices. 



Here we explore theoretically certain properties of cir- 
cular graphene quantum dots, defined via cutting the de- 
sired shape from a two-dimensional extended sheet. In 
particular, we regard investigations of graphene QDs as 
providing an opportunity for reexamination (and possi- 
bly experimental resolution) of remaining questions con- 
cerning the appropriateness of liquid-type vs. molecular- 
type trial functions for a finite number of 2D electrons. 
Indeed, it has been known for some time that mani- 
folds of degenerate midgap edge states exist in graphene 
nanostructures (such as graphene ribbons) when they ter- 
minate in a zig-zag boundary. 36 ' 37 In a recent paper 38 
it was noted that the single-particle edge states associ- 
ated with circular graphene dots with zig-zag boundary 
conditions, and in the absence of an applied magnetic 
field, display degeneracies and quantum numbers in close 
analogy with the manifold of single-particle states that 
form the familiar LLL in semiconductor heterostructures 
at high B. Furthermore, the numerical calculations of 
Ref. 38, covering rather limited ranges of electron num- 
bers (i.e., 2 < N < 5) and total angular momenta (i.e., 
N(N—1) < L < 60), suggested that the use of quantum- 
liquid-type trial functions in relation to the graphene LLL 
(gd-LLL) may be less promising than that of Wigncr- 
crystal-type ansatzes. 



In this paper, applying a methodology based on 
angular-momentum projection techniques that was intro- 
duced in Ref. 8, we derive analytic REM trial functions 
appropriate for the gd-LLL. By introducing a single vari- 
ational parameter, we demonstrate numerically (through 
systematic comparisons with EXD calculations) that the 
variational variant of the REM (referred to as vREM) 
substantially outperforms the Laughlin trial functions (as 
well as the ansatz of Ref. 38) for all values of fractional 
fillings within the expanded angular-momentum range 
N(N — l)/2 < L < 150, and for all of the following 
larger numbers of electrons N — 5, 6, 7 and 8. 



II. LOWEST LANDAU LEVEL FOR CIRCULAR 
GRAPHENE DOTS 

A. Single-particle edge states 

It is well known 39 that the low-energy bandstructure of 
graphene can be described by a linearized tight-binding 
Hamiltonian H . For a graphene dot with a circular sym- 
metry this linearized Hamiltonian is given 38 in polar co- 
ordinates by 

H=hv F ( H + £_y (i) 

where 

H - ( ^ (~ i9 r ~ 7 9 *) \ ( \ 

s ~\ (-id r + id*) ) ' {Z) 

where vf is the Fermi velocity, and s = ± specifies the 
degenerate in energy valleys for the two bands formed 
at the K and K 1 points. The index s can be considered 
as a "pseudospin" , which creates a fourfold degeneracy 
when the spin degree of freedom is also considered. The 
general solution of the Hamiltonian in Eq. (2) is a two 
component vector of the form 

( u?(r,4>) \ /ox 
{uf(r,4>))> (3) 

where A and B denote the two triangular sublattices of 
graphene. 

The usual volume solutions (which are zero on the 
graphene boundary but otherwise extend everywhere in- 
side the area enclosed by the graphene dot) have energy 
Ek = vpk, with uf and uf components that are ex- 
pressed via the Bessel functions. Here we are not inter- 
ested in such volume solutions. Instead we focus on the 
special edge states with zero energy E = 0. These E = 
states are eigenfunctions of H s under the assumption that 
the graphene boundary exhibits an uninterrupted zigzag 
edge; 36-38 an outline of their derivation from the Hamil- 
tonian H s is given in Appendix A. 

Henceforth we will only need to remember the precise 
form of the edge states, which is given by 

and 

( %~ ) = ( v^S^ ) ' (5) 

Namely one of the A and B components is everywhere 
zero (both on the boundary and inside the dot) and the 
two valleys + and — are decoupled even when the two- 
body Coulomb interaction is considered (which is the 
main focus of this paper; see below). As a result, in 
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the following, we will drop the sublattice and valley in- 
dices. We will also assume that the electrons are fully 
polarized. 

Since the single-particle angular momentum I > (to 
guarantee normalizability) , the manifold of such model 
edge states forms a set of degenerate states similar to 
the lowest Landau level (LLL), familiar from the case of 
2D semiconductor devices at very high magnetic fields 
B. We will call the manifold 40 of degenerate edge states 
with I > the graphene-dot lowest Landau level (gd- 
LLL). The main difference [apart from the normalization 
constant, see Eq. (3) in Ref. 8] between the two cases is 
that the single-particle states in the usual LLL exhibit an 
additional Gaussian multiplicative factor exp(— r 2 /4A|) 
where As = Tic/(eB) is the magnetic length. This Gaus- 
sian is missing from the expression for the edge states in 
Eqs. (4) and (5); instead one has tpi = for r > R. 

It is thus natural to investigate possible similarities 
related to fractional quantum Hall effect (FQHE) physics. 



B. Classes of variational many-body wave functions 

FQHE physics in the LLL has been extensively investi- 
gated for the case of 2D semiconductor quantum dots. 7 ' 11 
A main focus has been the underlying nature of the cor- 
related many-body states, i.e., 'liquid' (Laughlin, com- 
posite fermions) or molecular ('crystalline', REM). De- 
tailed comparisons of pair-correlations functions between 
JL/CF and REM states with EXD ones support the view 
that the molecular (localized electrons) picture in semi- 
conductor QDs provides the most appropriate descrip- 
tion. The emergence of a gd-LLL, as described above in 
graphene dots offers a further area for testing the ap- 
propriateness of liquid-type variational wave functions 
(JL/CF) versus those that describe REMs. 

First we will proceed with deriving a modified REM 
trial wave function that takes into consideration the dif- 
ferences between the single-particle states which span the 
usual LLL (zero-node 2D-harmonic-oscillator states) and 
the gd-LLL (edge states). 



III. DERIVATION OF VARIATIONAL REM 
TRIAL WAVE FUNCTIONS FOR GRAPHENE 
DOTS 

A. Intermediary parameter-free REM functions 

REM analytical wave functions in the LLL for elec- 
trons in two-dimensional semiconductor quantum dots 
were derived earlier in Ref. 8. The physics underlying 
such a derivation is based on the theory of symmetry 
breaking at the mean-field level and of subsequent sym- 
metry restoration via projection techniques. 11,41 In par- 
ticular, this approach consists of two steps: 

(I) At the first step, one constructs a Slater deter- 
minant 1 $ N (zi, . . . , zn) out of displaced single-particle 




FIG. 1: The displaced orbital u(z, Zo) (modulus square) rep- 
resenting a localized electron at the point Zo — 1 + Oi. The 
radius R of the dot serves as the unit of length. 



states u(zj, Z .j), j = 1, . . . , N, that represent the elec- 
trons localized at the positions Z j, with (omitting the 
particle indices) z = x + iy = re 1 ^ and Zo = Xo + iYq = 
Roe 1 ' 1 ' . Note that necessarily all electrons are localized 
radially on the edge of the graphene dot, so that i?o = R- 
Naturally, for the LLL case of semiconductor QDs, the 
localized u(z, Zo) single-particle states (referred to also 
as orbitals) were taken to be displaced Gaussians with 
appropriate Peierls phases due to the presence of a per- 
pendicular magnetic field [see Eq. (I) in Ref. 8]. In the 
case of electrons in graphene dots, however, the gd-LLL 
is spanned by edge-like orbitals (without a Gaussian fac- 
tor), i.e., 



1 + 1 z l 
^RZR 1 ' 



(6) 



and result the appropriate localized orbitals arc 

taken to have an exponential form 



u(z,Z ) = Gexp(z/Z ), 



(7) 



with G being the normalization constant (depending only 
on R). The fact that u(z,Zo) in Eq. (7) represents a 
localized electon is illustrated in Fig. 1. 

The localized orbital can be expanded in a series over 
the basis functions in Eq. (6) in the following way 



u(z,Z ) = J2 C i( Z °)Mz) 



Z=0 



with 



Ci(Zo) = G 



^R l+1 1 

nVT+iz ' 



(8) 



(9) 



When constructing the many-body Slater determinant 
\E ,JV [,z], one considers N orbitals u(zj,Z j) representing 
N electrons on a ring of radius R (the radius of the 
graphene dot) forming a regular polygon, i.e., 



Zo.j = Re 2m(1 - j ^ N , l<j<N. 



(10) 
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The single Slater determinant f^Jz] represents a static 
electron (or Wigncr) molecule (REM or RWM). Using 
Eq. (8), one finds the following expansion (within a pro- 
portionality constant): 



A'r 



E 



y/{h + 1) . . . (l N + 1) 

Mw 



h,...,l N =0 

xC h (Z .i) ■■■Ci N {Z 0iN )D(l u l 2 , l N ), (11) 



Ji 



the ele- 
with 



where D(h,l 2l ■ ■ ■ Jn) = det[z[\ z l 2 2 , . . . , z l £ 
mcnts of the determinant are the functions zs 
Zi, z l 2 2 , . . . , z 1 ^ being the diagonal elements. 

(II) Second step: The Slater determinant ^/^[z] breaks 
the rotational symmetry and thus it is not an eigenstate 
of the total angular momentum Tit = Ti^^Li h- How- 
ever, one can restore ' 11 the rotational symmetry by ap- 
plying onto ^ 



N 



z] the projection operator 

2tt 



V l = 



2Wo 



d'ye 



(12) 



where TiL are the eigenvalues of the total angular mo- 
mentum. 

When applied onto ^^[2], the projection operator Vl 
acts as a Kronecker delta: from the unrestricted sum in 
Eq. (11) it picks up only those terms having a given total 
angular momentum L (henceforth we drop the constant 
prefactor Ti when referring to angular momenta) . As a re- 
sult the projected wave function <j>^ = 
as (within a proportionality constant) 

JiH — Hn=L 



Vl^ is written 



N r 



Z = 



E 

h,...,h 



D(Ii,...,In) j(4>°h+-+^° N i N ) 
h\...l N \ 



(13) 



with </>° = 2tt(j- l)/N. 

We further observe that it is advantageous to rewrite 
Eq. (13) by restricting the summation to the ordered ar- 
rangements l\ < I2 < •■• < In, m which case we get 



Nr 



h+h^ Hn=L 

E 

0<h<l 2 <-<lN 



xdet[e^ Zl ,e^' 2 , 



D(k,...,l N ) 
h\...l N \ 

i4>° N l N ] 



(14) 



The second determinant in Eq. (14) can be shown 42 
to be equal (within a proportionality constant) to the 
following product of sine terms times a phase factor (in- 
dependent of the individual Z/s): 



in(N-l)L/N 



n 



sin 



h) 



(15) 



l<j<k<N 

Thus, the final result for the REM wave function is 
(within a proportionality constant): 



REM r 
N,L I 



E 

o<ii<; 2 <---</iv 
x Jf sin 

l<j<k<N 



D(l\,h, ■ • ■ , In) 



hU 2 \. 

7T 
N 



..l N \ 
Ik) 



(16) 



B. Introducing the variational parameter 

As described below, we found that the agreement be- 
tween the REM in graphene dots and the EXD solutions 
can be improved in a nontrivial way by introducing vari- 
ational parameters. In particular, we found that con- 
sideration of a single variational parameter a serves our 
purpose remarkably well. Specifically, one replaces the 
prefactor 



(17) 



h\l 2 \...l N \ 



in Eq. (16) by the following expression: 

[(Zl + l)(/2 + l)...(^ + l)] (1 - a)/2 



(h\l 2 \...l N l) 



no 



(18) 



We call the a-optimized wave functions the variational 
REM functions (denoted by vREM). When a = 1, the 
vREM coincides with the parameter-free REM expres- 
sion. We note that a single-parameter variational crystal- 
type wave function, but with a different dependence on 
the parameter a, has also been employed in Ref. 38. The 
present choice of variational parameter [see Eq. (18)] pro- 
duces substantially better results (see below). From a 
practical point of view, we note that the crystal-type 
wave function proposed in Ref. 38 does not contain a 
"less-than" ordered-arrangement restriction in the sum- 
mation indices Ii,...,In, and as a consequence it gener- 
ates an exponentially larger number of expansion terms, 
thus greatly inhibiting numerical evaluations for larger 
N and L. 



IV. EXACT DIAGONALIZATION AND 
TWO-BODY COULOMB MATRIX ELEMENTS 

For a circular graphene QD comprising N electrons 
in the gd-LLL, the many-body hamiltonian 7i comprises 
only the two-particle interelectron Coulomb repulsion, 
i.e., 



N N 



^-EE 



(19) 



i=l j>i 



where n is the dielectric constant and denotes the 
relative distance between the i and j electrons. 

The REM wave functions $^ E ^ 1 derived in the previ- 
ous section will be compared to the EXD ones 3>jvjP that 
are solutions of the exact diagonalization of the hamilto- 
nian (19) in the many-body Hilbert space spanned by the 
Slater determinants 



N\ 



ip h (z N ) ... ipi N (z N ) 



(20) 
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where the single particle functions (z) are given by the 
edge states of Eq. (6) and the index / counts the arrange- 
ments < l\ < h < ■ ■ ■ < In with l± + 12 + ■ ■ ■ + In = L. 
Namely, & n X l is written as 



(21) 



and the exact diagonalization of the many-body 
Schrodinger equation 



(22) 



yields the cocmccients Cl and the EXD eigenenergies 

pEXD 
^N,L ■ 

The matrix elements (5 , £|7i|\E , £) between the basis de- 
terminants [see Eq. (20)] are calculated using the Slater 
rules 43 and taking into account that, in the gd-LLL, 
the many-body hamiltonian has contributions from the 
Coulomb interaction only, i.e., 



%<3 



(23) 



Naturally, one also needs the two-body matrix elements 
of the Coulomb interaction in the basis formed out of the 
single-particle edge states. These matrix elements are 
given through appropriate analytic expressions. Indeed 
by defining 



M(m,n,k) = 
dzi / dz 2 ^ m+k (zi)^ k {z2)- 



z\ - z 2 \ 



(24) 



finds 



M(m, n, k) = 



1 2ttC0F r (k 



R(2m + 2n + 3) T(k + l) 



(25) 



r(n+ 1) /1/2, fc + 1/2, re + 1 
r(n + 2) 3 2 V fc + 1, n + 2 



1 



T(m + fc + l) 
T(m + k + 2) 



3F2 



1/2, A; + 1/2, m + k+1 
k + 1, m + k + 2 



where 3F2 is the generalized hypergeometric function 44 
at the point x = 1, T is the Gamma function, and 

c = y/(™, + k + 1) (n - k + 1) (m + 1) (re + f) 



V. NUMERICAL RESULTS 

A. EXD total energies 

In Fig. 2 we display systematic EXD total energies 
in the range of TV = 5 to TV = 8 edge electrons as a 
function of the total angular momenta L (in the large 
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FIG. 2: Exact diagonalization ground-state energies in the 
graphene-dot LLL for N = 5, 6, 7, 8 electrons, as a func- 
tion of the total angular momentum L. Observe the appear- 
ance of cusp states of enhanced stability at the magic angular 
momenta L m = N(N - l)/2 + kN , k = 0, 1, 2, . . ., a fact 
that indicates formation of Wigner molecules having a sin- 
gle polygonal-ring configuration (0, N). Energies in units of 
e 2 /(K,R), with k being the graphene dielectric constant and 
R the radius of the quantum dot. For L — > 00, the ground- 
state energies approach asymptotically the classical electro- 
static energy [see Eq. (27)]. 



range < L < 150). This large L range and the con- 
sideration of N > 5 electrons were not reached in an- 
other recent publication; 38 they are, however, essential 
for unequivocally establishing the proper similarities and 
differences with the high-magnetic-field physics of semi- 
conductor QDs. 

For fully polarized spins considered here, the minimum 
total angular momentum is Lq = N(N — l)/2, in analogy 
with the case of semiconductor QDs. 11 Furthermore, in 
analogy again with the case of semiconductor QDs, the 
total energies decrease on the average as L increases. On 
top of this average trend, one observes prominent oscil- 
lations of period N. These oscillations reveal that the 
states with L = Lq + kN, k = 0, 1, . . ., are energetically 
the most stable in their immediate neighborhood. Bor- 
rowing the terminology from the literature 11,45 of semi- 
conductor QDs, we refer to these states in graphene 
QDs as cusp states, and the corresponding total an- 
gular momenta (i.e., L = Lq + kN) as magic angular 
momenta. It is well known that cusp states develop 
to fractional-quantum-Hall-effect (FQHE) states in the 
thermodynamic limit (TV — > 00), with the corresponding 
fractional filling factor being v = Lq/L. 

Following a similar analysis 11,45 with the case of semi- 
conductor QDs, one can conclude that the appearance of 
the oscillatory period N in the total energies (associated 
with the cusp states) is a reflection of formation of (0, N)- 
type Wigner molecules, with all the electrons localized on 
a single ring (of radius R) at the apices of a regular TV- 
polygon. There is a major difference, however, between 
the present system and the semiconductor quantum dot 
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FIG. 3: Conditional probability distribution [see Eq. (28)] as- 
sociated with the EXD ground state in the gd-LLL of N = 7 
electrons and for L = 63 (corresponding to fractional fill- 
ing v = 1/3). One clearly observes 6 humps in agreement 
with formation of a (0, 7) Wigner molecule, and in contrast 
to the liquid-like Laughlin physical picture. The fixed (obser- 
vation) point is denoted by a solid dot. Lengths in units of 
the graphene-dot radius R. Vertical axis in arbitrary units. 



case. That is, in semiconductor QDs, more than one iso- 
mers may form with concentric multiring arrangements 
occurring for N > 5 electrons in the dot; such arrange- 



ments are denoted as (n-i, n,2, . 



. n, 



,) (see Ref. 11), where 



n r , r = 1,2, ... ,q are the number of localized electrons 
on each ring; X)r=i Ur = ^ • ^ n contrast, in the case of 
graphene dots only the one-ring (0, N) molecular configu- 
ration arises (with no electron residing at the geometrical 
center of the graphene QD). 

For L — > oo, the EXD energies in Fig. 2 approach the 
limiting value corresponding to the classical electrostatic 
energy of N point-like electrons in a (0,N) configuration 
with radius R, i.e., 



E cl {N) 



4kR 



NS N , 



(27) 



with S N = Ef =2 (sin[(j - l)^/N}y\ 



B. EXD densities and pair correlations 

The EXD eigenfunctions conserve the total angular 
momentum and the corresponding electron densities are 
circularly symmetric. This property "conceals" the pres- 
ence of the Wigner molecule. The crystalline structure, 
however, is present in the intrinsic frame of reference of 
the electron molecule, and it can be revealed through 
the use of the fully anisotropic pair correlation function 
P(z,zq), defined as 



P(z, z ) 



(* E N X F\ 



6(z - Zi)5(z 



Zj )\$l x E). (28) 



P(z, Z q) is often referred to as the conditional proba- 
bility distribution (CPD), since it is proportional to the 



FIG. 4: Relative error (per electron) of the vREM ground- 
state energies as a function of the total angular momentum 
L. 



probability of finding an electron at z under the condition 
that another one is situated at the point zq (the socalled 
fixed, or observation, point). 

In Fig. 3 we display the conditional probability distri- 
bution for the case of N = 7 electrons and the magic total 
angular momentum L = 63 (L — Lq = 7k, with k = 6), 
which corresponds to the celebrated v = 1/3 fractional 
filling. One clearly observes six humps (arranged in a 
single-ring configuration) associated with formation of a 
(0,7) rotating Wigner molecule. (As is well known from 
the literature of semiconductor quantum dots, 11 the lo- 
calized electron at the fixed point does not contribute any 
hump in the CPDs.) Similar CPDs are found for other 
values of N. 



C. Comparison between vREM and EXD wave 
functions 

We turn now to comparisons between the EXD wave 
functions and the vREM ones. We first observe that the 
REM and vREM functions [see Sect. Ill] correspond to 
the magic angular momenta L = Lo + kN, k = 0, 1, 2, . . ., 
since all the sine-product coefficients in the expansion 
(16) are identically zero for L / Lq + kN. In this sec- 
tion, we will show that the vREM functions represent 
a high-quality approximation to the EXD eigenfunctions 
by investigating wave function overlaps and relative er- 
rors between the total energies obtained by the two meth- 
ods; the relative errors are defined as A v rem(-^, L) — 

jvREM _ £'EXD^ / pEXD 



(E 



N.L 



N,L 



)/E 



N.L 



We start by displaying in Fig. 4 the relative error 
A vREM (N,L) of the vREM total energies. The vREM 
offers an excellent approximation, since the maximum 
relative error is less than 0.045%. For all sizes exam- 
ined, the maximum relative error occurs about v = 1/3 
(see Section V A) , and subsequently it decreases as L in- 
creases, approaching zero as L — > oo. 

In Fig. 5, we display the overlaps <S v rem = 



7 




($EXD| $ vREM) between the vREM functions and the 
EXD solutions. These overlaps are larger than 0.985 for 
N = 5 and larger than 0.95 for N — 8 and they tend to 
slowly approach unity as L increases. 

In Fig. 6, we display the values of the variational pa- 
rameter a that optimize the vREM trial functions for a 
given number of electrons N as a function of L. These 
values are significantly different from unity (which corre- 
sponds to the parameter-free REM). In fact the optimal 
a values are smaller than 0.4, and they slowly decrease to 
about 0.16 for L = 150, for all the values 5 < N < 8. We 
stress that optimization of a is essential for achieving 
a high quality reproduction of the EXD ground states. 
Without the additional optimization (i.e., taking only 
the value a = 1) the behavior of the overlaps is un- 
satisfactory, since they tend to diminish as L increases 
(see Appendix B). The degradation of the overlaps of the 
parameter-free REM (a = 1) as L increases in the case 
of the graphene quantum dot contrasts with the oppo- 
site behavior of the overlaps of the parameter-free REM 



in the case of semiconductor quantum dots. 8,9,11 This 
difference is attributed to the absence of translational in- 
variance for the electrons in the graphene quantum dot, 
which leads to differences in the organization of the EXD 
excitation spectra. 



D. EXD versus Laughlin wave functions 

It is interesting to compare the accuracy with which 
the vREM wave functions approximate the EXD ones 
with that of the Laughlin trial functions. The Laughlin 
wave functions are restricted to the socalled main (odd) 
fractions v = l/(2m + 1) and have played an impor- 
tant role in the FQHE literature of the extended two- 
dimensional electron gas in semiconductor heterostruc- 
tures. Their form is 

®NT n M= II (^-^) 2m+1 : ( 29 ) 
l<i<j<N 
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FIG. 6: The values of the variational parameter a [see Eq. 
(18)] that optimize the vREM trial functions for a given N as 
a function of L. 



where the Gaussian factors are missing [see Sect. II A] 
due to the differences in the single-particle states be- 
tween semiconductor and graphene quantum dots; m = 
0,1,2,.... 

In Fig. (7), we display the relative error, 
A Laughlin (iV,i) = - EfpyEf™, of 

the Laughlin total energies with respect to the ground- 
state EXD ones as a function of increasing total angular 
momentum L. The Laughlin relative errors are substan- 
tially larger (on the average by a factor of 5) than the 
vREM ones [see Fig. 4]; this is the case even for the 
celebrated v = 1/3 fractional filling. 

In addition, the Laughlin overlaps <SLaughiin = 
($EXD| $ Lau g hii„^ plotted in Fig 8 ) exnibit an unsat- 
isfactory performance compared to that of the vREM 
overlaps, that is: (i) they become steadily smaller as the 
angular momentum increases, and (ii) even for v = 1/3, 
they are smaller than the corresponding vREM overlaps 
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in all instances studied here, i.e., TV — 5 — 8 electrons in 
the graphene dot. 

We conclude that the Laughlin functions fail to capture 
the case of the gd-LLL, while the vREM functions offer 
an appropriate approximation for graphene QDs. 



E. EXD versus composite-fermion wave functions 

It is also interesting to compare the accuracy with 
which the vREM wave functions approximate the EXD 
ones with that of the composite-fermion trial functions, 
which are more general than the Laughlin functions. 
Along with the Laughlin functions, they CF trial func- 
tions have played a significant role in the FQHE literature 
of the extended two-dimensional electron gas in semicon- 
ductor heterostructures. Their form 5,7 in the disc ge- 
ometry (case of 2D QDs studied here) are given by the 
expression, 

^ F (N)=V LLL J] { Zi - Zj f m V™, (30) 

l<i<j<JV 

where z = x+iy and ^ l ^ t M is the Slater determinant of TV 
non-interacting electrons of total angular momentum L* ; 
it is constructed according to the Independent Particle 
Model (IPM) from the Darwin- Fock 46 orbitals ippj{z), 
where p and I are the number of nodes and the angular 
momentum, respectively [for the values of p and I in the 
nth Landau level in high B, see Appendix F]. 

The single-particle electronic orbitals in the Slater de- 
terminant ^^, M are not restricted to the lowest Landau 
level. As a result, it is necessary to apply a projection 
operator Vlll to guarantee that the CF wave function 
lies in the LLL, as appropriate for B — > oo. We carry the 
T^lll projection according to section 4.3 of Ref. 47. Af- 
ter obtaining the projected CF function in the LLL, the 
corresponding trial function in the gd-LLL is constructed 
by simply replacing ip^J(z) by the ipi(z) in Eq. (6). 
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FIG. 8: Overlaps of the Laughlin trial states with the EXD 
ones as a function of total angular momentum L. 



TABLE I: Total CF and EXD energies [per particle, in units 
of e 2 /(fti?)] for TV = 6 electrons in a graphene dot. The 
corresponding relative errors, Acf/N, and overlaps, <Scf, are 
also listed. For the determination of the auxiliary angular 
momenta L* , see Appendix F. 



L(L") 


E EXD /N 


E CF /N 


10" 4 Acf/N 


ScF 


21(-9) 


2.31357 


2.3548 


29.67 


0.887 


51 (-9) 


1.99693 


2.0385 


34.67 


0.793 


30(0) 


2.24863 


2.3473 


73.17 


0.369 


60(0) 


1.99452 


2.0520 


48.00 


0.507 


35(5) 


2.15628 


2.3022 


112.8 


0.356 


65(5) 


1.97413 


2.0410 


56.50 


0.451 


39(9) 


2.06465 


2.0952 


24.67 


0.892 


69(9) 


1.94477 


1.9713 


22.67 


0.754 



Since the CF wave function is an homogeneous poly- 
nomial in the electronic positions Zj's, its angular mo- 
mentum L is related to the non-interacting total angular 
momentum L* as follows, 

L = L* +mN(N- 1) = 2mL . (31) 

Here we will consider the mean-field version of the 
composite-fermion theory, according to which the Slater 
determinants \E'^ P , M are the socalled compact states (see 
Appendix F for details; the corresponding values of L* 
are listed in Table II). We note that recently several ex- 
tensions of the CF theory have been formulated 48 that 
account for residual-interaction effects among the indi- 
vidual $1 F (N) composite-fermion states. Consideration 
of such residual-interaction effects is beyond the scope of 
the present paper. 

In Table I, we compare total CF and EXD ener- 
gies (per particle) for TV = 6 electrons in a graphene 
dot. We also display the corresponding relative errors 
A CF (TV, L)/N = (E%\ - Ef™)/{NEff) and overlaps 

In addition to the L = 15 + 6k, k = 0, 1, 2, . . . magic an- 
gular momenta for TV = 6 found from EXD calculations 
(see Fig. 2), the compact-state CF theory mistakenly 
predicts the existence of magic angular momenta with 
L = 15 + 5fc, k = 0, 1, 2, . . ., e.g., for L = 30, 35, 60, 65. 
Furthermore, even for the states with L = 15 + 6fc, e.g., 
L = 21, 39, 51, 69 (Table I), the quantitative performance 
of the compact CF functions (concerning relative errors 
and overlaps) is rather weak compared to that of the 
vREM: the CF relative errors are larger roughly by a 
factor of 10, while the CF overlaps are systematically 
smaller (< 0.9) than the vREM ones (> 0.97) (see Figs. 
2 and 5, and Table I). 

As was the case with the Laughlin functions, we con- 
clude that the compact CF functions are also at a disad- 
vantage compared to the vREM concerning the descrip- 
tion of strongly correlated states in the gd-LLL. 
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FIG. 9: Relative errors (per electron) of the energies of the 
WC-ansatz [Eq. (32)] as a function of total angular momen- 
tum L. 



Comparison with the Wigner-crystal ansatz of 
Ref. 38 



Here we compare the vREM total energies with those 
associated with the Wigner-crystal trial function of Ref. 
38, given by 



<I> WC M - 



h+h-l \-In=L 



exp 



'£ 



2'Knl n 

N 



x]l(l n + l) w+1/2 D(h,l 2 ,...,l N ), (32) 



where w is a variational parameter. 

In Fig. 9, we display the relative errors Awc(-W, £) = 
{Ejj c L -E^)/E^ of the WC-ansatz energies relative 
to the EXE) ones. From a comparison of these results 
with those displayed in Fig. 4, we conclude that the rel- 
ative errors of the WC ansatz are on the average at least 
twice as large as those corresponding to the vREM, re- 
flecting the superior description of the gd-LLL provided 
by the latter function. 



VI. PINNED ELECTRON MOLECULES 

The zig-zag geometry (on which the boundary condi- 
tions are applied; see Section II A) does not allow forma- 
tion of a continuous circular edge without some structural 
or chemical modification of the graphene hexagonal lat- 
tice structure. Without such modification, regions along 
the circular edge satisfying a zig-zag condition must nec- 
essarily be disrupted by a number of discrete points asso- 
ciated with arm-chaired conditions. 49 It has been found 
that the edge states are robust in this case, 37,49 and as 
a result these discrete set of disruptions act as effective 
impurities that modify the many-body hamiltonian in 
Eq. (19). The presence of such impurity terms in the 
many-body hamiltonian will mix the good-total-angular- 
momentum REM states, resulting in the formation of 



FIG. 10: Electron density of a pinned molecule for N = 7 
electrons formed from the linear superposition of two REM 
states with L = 35 and L — 42. Lengths in units of the 
graphene dot radius R. Electron density in units of R~ 2 . 



pinned electron molecules (PEMs). In contrast to the 
REMs (whose electron density is uniform along the az- 
imuthal direction, that is, not showing any azimuthal 
density modulation), the electron density of a pinned 
electron molecule is expected not to have circular sym- 
metry; it will exhibit angular density oscillations, and the 
number of humps will equal the number of electrons N. 

We demonstrate this property of a PEM for two partic- 
ular cases displayed in Figs. 10 and 11. Fig. 10 displays 
for N = 7 the electron density for the linear superposi- 
tion of two REM states with L = 35 and L = 42, while 
Fig. 11 displays for N = 8 the electron density for the 
linear superposition of two REM states with L — 44 and 
L = 52. In both cases the expected angular modulation 
is clearly well formed with seven humps in the former 
and eight humps in the latter case. 



ED N=8 (L=44)+(L=52) 




FIG. 11: Electron density of a pinned molecule for N = 8 
electrons formed from the linear superposition of two REM 
states with L = 44 and L — 52. Lengths in units of the 
graphene dot radius R. Electron density in units of R~ 2 . 
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VII. SUMMARY 

The manifold of degenerate midgap (zero-energy) edge 
states in circular graphene quantum dots with zig-zag 
boundaries resembles, under free- field conditions, 50 the 
celebrated lowest Landau level, familiar from the case of 
semiconductor heterostructures in high magnetic fields. 
The effect of e — e interactions in this graphcnc-LLL 
were systematically investigated and were found to gener- 
ate many-body strongly correlated behavior that exhibits 
many similarities with the fractional quantum Hall effect. 

Numerical exact-diagonalization studies were pre- 
sented for 5 < N < 8 fully spin-polarized electrons and 
for total angular momenta in the range of N(N — l)/2 < 
L < 150. Moreover, we presented a derivation of a 
rotating-electron-molecule type wave function based on 
the methodology introduced earlier 8 in the context of 
the FQHE in two-dimensional semiconductor quantum 
dots. The EXD wave functions were compared with the 
derived rotating-electron-molecule and other suggested 
FQHE trial functions, like the Laughlin function and the 
Wigner-crystal ansatz of Ref. 38. It was found that a 
variational extension of the REM offers a better descrip- 
tion for all fractional fillings compared with that of the 
Laughlin and Wigner-crystal ansatz functions (including 
total energies and overlaps). The success of the REM 
function reflects the importance of strong azimuthal lo- 
calization of the edge electrons in graphene quantum 
dots. 

The variational REM functions were derived through 
the use of a two-step method: (i) first a mean-field-type 
single Slater determinant constructed out of N localized 
electron orbitals (that break circular symmetry) was con- 
sidered; this determinant describes the finite analog of 
a classical static Wigner-crystal, and (ii) a multideter- 
minantal wave function was generated through the sub- 
sequent application of projection techniques that intro- 
duced azimuthal fluctuations and restored the circular 
symmetry and good total angular momenta. 

In contrast with the multiring arrangements of elec- 
trons in circular semiconductor quantum dots, we found 
that the graphene REMs exhibited in all instances a sin- 
gle (0, N) polygonal-ring molecular structure. Disrup- 
tions in the zig-zag boundary condition along the circu- 
lar edge behave effectively as crystal-field effects that pin 
the electron molecule, yielding single-particle densities 
with broken rotational symmetry that portray directly 
the azimuthal localization of the edge electrons. 
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FIG. 12: Overlaps of the parameter-free trial REM states [i.e., 
for a = 1, see Eq. (16)] with the EXD ones, as a function of 
total angular momentum L. 



APPENDIX A: MORE ON EDGE STATES 

The general solution of the eigenvalue equation corre- 
sponding to the linearized tight-binding Hamiltonian (2) 
is of the form 



Bs 



WW r ) e i[i+(l-s)/2]<?S 
x Bs^ e i[l+(l+s)/2]4, 



(Al) 



where s = ±; obviously s — ±1 when occurring in a 
phase. As a result, the matrix eigenvalue problem is 
equivalent to the following set of equations involving the 
vector components: 



-ihvpd r xf + — ihvF{l + 1) 
—ihvFd r xf + + ITivfI 



Xi 



B + 



A+ 



X, 



= E X f + , (A2) 



where we considered only the case for s = + (the s = — 
case can be treated in a similar way). 

We are interested in solutions with E = (the socalled 
midgap solutions), in which case the set of equations (A2) 
reduces to 



d rX r 



A+ 



-Xi 



r Xl 



B+ 







= 0. 



(A3) 



The solutions of these equations are 
*f+(r) = Xf + (R)Q 



X? + (r) 



X, 



A + 



(A4) 



The boundary condition is that of a zigzag graphene 
edge that ends always on a site of the same lattice, i.e., 
the condition 



(A5) 
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FIG. 13: Overlaps of the variational REM states with the 
EXD first-excited states for selected values of the variational 
parameter a [see Eq. (18)] as a function of total angular mo- 
mentum L. The drops toward zero determine the optimal a's 
for given L's. For a = 1 (top curve), no such drop to zero 
occurs. 

forces the B+ component to vanish everywhere on the B 
sublattice, yielding the final form 

(^1) = («*>^). m 

The normalization constant xt + {R) ' IS easily calcu- 
lated and was given in Eq. (4). 



J2i Cl^lI*]' where the basis wave functions \E'£ are the 
Slater determinants defined by Eq. (20). 
Specifically one has 

N 

p{z) = ($ L \J2S(z-Zk)\$L) 

k=l 

/ ./ fc=i 

N 

/ fe=i 

AT 

= Ei^i'E^w^^) 

I k=l 

/ fe=i 

where the edge states ^'s are given by Eq. (6) and l J k 
denotes the single-particle angular momenta associated 
with the Slater determinant Vf^; naturally L = J2k=i 'fe- 
The single-particle density operator connects in principle 
Slater determinants that differ at most in one orbital. 
However, in the LLL, the conservation of the total an- 
gular momentum implies that there is no pair of Slater 
determinants in the linear superposition of that differ 
precisely by one single orbital; thus one sets J = I when 
deriving Eq. (CI). 



APPENDIX B: THE PARAMETER-FREE REM 

(a = l) 

As mentioned in Section VC, the overlaps 6>rem = 
(^iv^l^ATL 1 ) between the parameter- free REM waves 
function in the gd-LLL and the EXD ones behave in an 
unsatisfactory way, i.e., they decrease as L increases. The 
precise behavior of Srem is displayed in Fig. 12, and 
it contrasts with that of the variational REM, <S v rem, 
displayed in Fig. 6. 

The degradation of the Srem reflects the fact that pro- 
gressively the overlap of the REM wave function with the 
excited EXD states increases with increasing L. On the 
other hand, the optimized values of a correspond to vari- 
ational REM trial functions that have practically zero 
overlap with these excited EXD ones. We have found 
that such optimal a values can be found for all studied 
values of N and L. This is illustrated in Fig. 13, where 
the overlaps of the vREM with the first excited EXD 
state dip toward zero at the optimal a values. 

APPENDIX C: SINGLE-PARTICLE DENSITY 

We give here the expression for calculating the single- 
particle density p{z) for a single many body state $£[2] = 



APPENDIX D: SINGLE PARTICLE DENSITY 
FOR A SUPERPOSITION OF TWO WAVE 
FUNCTIONS 

In Section VI we discussed how the disruptions in 
the zigzag boundary conditions create crystal-field effects 
that pin the rotating electron molecule. The effect of this 
pinning is described through the linear superposition of 
two many-body wave functions (EXD and/or REM) with 
magic good total angular momenta L and M; namely, 
through a wave function $ PIN such that 

$ PIN = -L (* L ± $ M ) , (Dl) 

where we have dropped the subscript N and superscripts 
'EXD' or 'REM' from the $z,(M)'s on the r.h.s. 

The $l(m)'s are known through their expansions over 
Slater determinants [see Section IV], i.e, 

$ L =^C£*£ and <&m = E C m*m, (D2) 
/ J 

and the Slater determinants ^ and "^ J M are built out of 
single particle states having individual angular momenta 
l{ and rrij such that J2k=i = L anQl SfcLi m k = 
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Using the operator p defined by the first line of Eq. 
(CI), the single-particle density is given by 

p plN (z) = ($ PIN |p|$ PIN ) 

= ^(($l|/S|$l) + ($m|/3|$m)± 

($l|/3|$m>±<*m|/3|$l». (D3) 



The diagonal terms ($l|/3|$l) and (3>m|/3|$m) are given 
by Eq. (CI). Since p is a one-body operator, the cross 
terms connect Slater determinants that differ precisely by 
one of orbital; 43 we denote by I 1 and m J q the correspond- 
ing pair of indices. By applying the Slater rules described 
in Ref. 43 (including bringing the two determinants into 
"maximum coincidence"), one finds: 



/ k=l J k=l 

^ 7 * T , ^ + 1)K + 1) ,r\"^ J 
± Y^Cl C J M a(L,I,p;M,J,q)V — 2 (-) 

± E C m C J L a(M,I,q;L,J,p)V _ (1) 



(D4) 



where o~{L, I,p; M, J,q) = ±1 depending on the even 
or odd number of exchanges of two rows (or columns) 
needed to bring the two determinants into maximum co- 
incidence. 



APPENDIX E: TWO-PARTICLE CONDITIONAL 
PROBABILITY DISTRIBUTION 

For the conditional probability density [see Eq. (28)], 
one has 

P(Z,Z ) = ££CM<*£|T|*£> ( E1 ) 
/ J 

where the operator T is symmetrized; 
T=X>- *0*(*o - Zj) + 5(z - Zj )S(z - sn). (E2) 

i<j 

The matrix elements of T between the two Slater de- 
terminants and ^ J L are calculated according to the 
Slater rules for a two-body operator. 43 

APPENDIX F: MORE ON COMPOSITE 
FERMIONS 

There is no reason to a priori restrict the Slater de- 
terminants ^> 1P * M to a certain form. 51 Following Ref. 
51, we will restrict the non-interacting L* to the range 
— Lq < L* < Lq, and we will assume that the Slater de- 
terminants , f^ p « M are the so-called compact ones. Let N n 
denote the number of electrons in the nth Landau Level 
(LL) with X)«=o N n — N; t is the index of the highest 



occupied LL and all the lower LL's with n < t are as- 
sumed to be occupied. The compact determinants are 
defined as those in which the N n electrons occupy con- 
tiguously the single-particle orbitals [if>® J(z)] of each nth 
LL \p + (\l\ — l)/2 = n] with the lowest angular momenta, 
I = — n, — n + 1, — n + N n — 1. The compact Slater 
determinants are usually denoted as [No, N\,..., N t ], and 
the corresponding total angular momenta are given by 
L* = (l/2)^l =0 N s (N s -2 S -l). 

For the CF theory, the magic angular momenta can be 
determined by Eq. (31), if one knows the non-interacting 
L*'s. For N = 6, the CF magic L's in any interval 
l/(2m - 1) > v > I /{2m + 1) [15(2m — 1) < L < 
15(2m + 1)], m — 1,2,3,4,..., can be found by adding 
2mLo — 30m units of angular momentum to each of 
the L*'s. To obtain the non-interacting L*'s, one needs 
first to construct 51 the compact Slater determinants. 
The compact determinants and the corresponding non- 
interacting L*'s are listed in Table II. 

There are nine different values of L*'s, and thus the 
CF theory for N = 6 predicts that there are always nine 
magic numbers in any interval 15(2m— 1) < L < 15(2to+ 
1) between two consecutive JL angular momenta 15(2m— 
1) and 15(2m + 1), m = 1,2,3,..., For example, using 
Table II and Eq. (31), the CF magic numbers for TV = 6 
in the interval 15<L<45(m=l) are found to be the 
following nine: 

15, 21, 25, 27, 30, 33, 35, 39. (Fl) 

In the interval 45 < L < 75 (to = 2), the CF magic 
numbers are: 

45, 51, 55, 57, 60, 63, 65, 69. (F2) 
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TABLE II: Compact non-interacting Slater determinants and 
associated angular momenta L* for N = 6 electrons accord- 
ing to the CF presciption. Both L* — —3 and L* = 3 are 
associated with two compact states each, the one with lowest 



energy being the preferred one. 


Compact state [No, Ni, N t ] 


V 


[1,1,1,1,1,1] 


-15 


[2,1,1,1,1] 


-9 


[2,2,1,1] 


-5 


[3,1,1,1] 


-3 


[2,2,2] 


-3 


[3,2,1] 





[4,1,1] 


3 


[3,3] 


3 


[4,2] 


5 


[5,1] 


9 


[6] 
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